function z1_fit = fit_std_dz_2state(parms,std_dz_target)

    [V,D] = eig(parms.StateTransitionProbs');
    [~,idx] = min(abs(1 - diag(D)));
    prob_z = V(:,idx)/sum(V(:,idx));
    
    opts = optimset('display','off');
    x_fit = fsolve(@fitfun,log(0.05),opts);
    z1_fit = -exp(x_fit);

    function z = fitfun(x)
        
        z = zeros(2,1);
        z(1) = -exp(x);
        
        condE_znext = parms.StateTransitionProbs*z;
        condE_dz = condE_znext - z;
        condV_dz = parms.StateTransitionProbs*(z.^2) - condE_znext.^2;
        var_dz = sum(prob_z.*condV_dz) + sum(prob_z.*(condE_dz.^2));
        
        z = sqrt(var_dz) - std_dz_target;
        
    end

end